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1. Introduction 

The billiard problem has played a vital role in the study of the manifestations of classical 
chaos in linear wave systems ("wave chaos") including microwave, optical and acoustic 
cavities and waveguides [UEllHlllllHlEllTllH], and various single-particle quantum 
systems [HI HOl [TT]. Even in strongly-interacting, non-linear systems the knowledge 
of the linear spectrum and eigenfunctions is paramount to infer complex observables 
[T2I IT31 [TH US, [Ml [IT]. In the semiclassical limit, or at high wavenumbers {k = 2it/X), 
the numerical solution of the Laplace eigenvalue problem becomes computationally 
challenging. Finite difference schemes [18] become impractical and Green's function 
matching methods suffer from the unfeasibility of a root search. 

The typical Green's function matching method (various implementations of which 
includes the Method of Particular Solutions (MPS) and boundary integral methods 
(BIM)) to solve the Laplace eigenproblem consists of finding the zeros of the secular 
determinant over a given wavenumber range. In practice, this is accomplished through 
the singular value decomposition (SVD) and scanning for the minima of the smallest 
singular values [20]. This requires typically of the order of (kR)^ matrix operations per 
mode (where R is the typical size of the system). Naturally, this procedure becomes 
progressively more expensive for higher lying eigenvalues. Missing eigenvalues are a more 
important problem in practice. At larger wavenumbers, when the spectrum becomes 
progressively denser, it's a serious problem to differentiate and separate the minima of 
the lowest singular valuesEl 

In this paper, we propose a fast and efficient method based on a Fredholm 
formulation of the billiard problem, to compute the spectrum and the corresponding 
eigenfunctions of the Laplace operator over a two-dimensional domain D. This method 
is closely related to the scattering quantization method (SQM) [211 1221 [23] as it relies on 
a similar acceleration technique of replacing the search for singular values of a matrix 
by an auxiliary eigenvalue problem. In contrast to SQM which expands the Laplace 
eigenfunctions in terms of a set of basis functions of the Laplace operator in the domain 
D, the expansion here contains the fundamental solutions of the Laplace operator. 
This has two important advantages which makes its exposition worthwhile. First, the 
proposed Fredholm formulation is known to be uniformly convergent [2l] while the SQM 
is known to be convergent only in so far as the Rayleigh hypothesis holds [25]. Second, 
Fredholm formulations via BIM are amenable to semiclassical quantization techniques 
through the transfer operator technique. Consequently, the behavior of the Laplace 
operator for various domain geometries in the semiclassical limit can be directly related 
to the invariants of underlying classical motion in that domain [261 [23, [28] . 

We would like to remark that the method outlined here provides a similar gain in 
speed and robustness with respect to the scaling method of Vergini and Saraceno [29l[30]. 
A recent boundary integral formulation of the scaling method has been carried out 

* Nearly degenerate levels can in practice be differentiated within the SVD scheme by looking at several 
of the smallest singular values [20] . 
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in [21] • To the authors' best knowledge the relation between scattering quantization 
methods and scaling methods is still an open question. 

2. Review of the BIM formulation 

Let us briefly review the BIM formulation of the billiard problem that we are addressing. 
Consider a two-dimensional Euclidean domain D bounded by a smooth boundary curve 
dD. Within this domain, let {'ip^{r)} be the set of eigenfunctions of the Laplace operator 
with eigenvalues kf^, 

-V'Mr) = klMr)- (1) 

We assume that ip^{r) satisfies Dirichlet boundary conditions ipfj,\dD = 0. In the context 
of the Schrodinger equation, E^j, = kf^ are the discrete energy levels of a particle in a box 
defined by dD. 

Associated to the differential operator in the above equation ([T]) is the Green's 
function 

{V^ + e)G{r,r';k) = 6{r-r'). (2) 

Regardless of the boundary conditions on the Green's function one can reformulate the 
billiard problem ([1]), through a completely standard procedure, in terms of a Fredholm 
integral equation of the second kind 

which has solutions only for discrete values k = k^. In the above equation, s is the 
arc length along the boundary, 8/ On = n{s) ■ V, d/dn' = n{s') ■ V, and n(s) is the 
outwards pointing unit normal of the boundary at location s. 




Figure 1. Schematics sliowing tlie variables used in tlie definition of tfie BIM kernel 
in equation ([5]). 

Therefore, the problem in the two-dimensional domain is reduced to a problem on 
the boundary. This reduction is physically very appealing as in the semiclassical limit 
the geodesic flow is uniquely represented as a discrete map on the boundary. Of course, 
the reduction in dimensionality has certain consequences. Whereas in the standard 
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treatment of domain problems through finite element methods one solves for the whole 
spectrum up to a maximal wavenumber k, boundary formulations provide a narrow 
spectral range around a reference wavenumber k. 

The standard BIM formulation employs the free-space outgoing Green's 
function [20l [32l [33] 

Go(r,r';fc) = -jH±(fc|r-r'|). (4) 

Here H^ (z) are the first and second kind Hankel functions of order zero. Let us rewrite 
the Fredholm problem ([3]) in an operator notation 

Ku = u, (5) 

where u{s) = -§^i^{f{s)) and the kernel, using the free Green's function (jlj) becomes 

= __cos^(.,/)H+(A:|r(.)-r(.')|). (6) 

Here, cos 6^(5, s') = n(s) ■ (r(s) — r(s'))/|r(s) — r(s')|, i.e. 9{s, s') is the angle between the 
normal at s and the cord connecting s and s' (see figured]). Consequently, K (referred 
to as K{k) in alternative notation) is clearly not a symmetric operator. Note that the 
diagonal elements are finite and given by 

lim K{s,s';k) = ^K{s), (7) 

where k{s) is the curvature at s. Hence the condition of quantization is 

det(l-ir(A;)) = 0. (8) 

The standard numerical procedure to extract the zeros of this secular determinant in 
the context of billiards is outlined in 



3. Scattering quantization approach to BIM 

In contrast to the standard procedure outlined in the last section, we shall reformulate 
the problem by considering the solution of the eigenvalue problem 

K{k)u = \u. (9) 

This eigenvalue problem provides us with a set of eigenvalues and eigenf unctions, 
{A''*^(/c), M^*-*(fc; s)} parametrically dependent on the continuous variable k. The structure 
of the operator K{k) is interesting. It can be shown via stationary phase integration 
that in the semiclassical limit, KK'^ is asymptotically diagonal i.e. while the off-diagonal 
elements are 0{-\/k), the diagonal elements are 0{k). The form of the diagonal elements 
is given by 
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Figure 2. (a) Distribution of the magnitude of the eigenvalues of K(}i) for kR = 20 
(red), 100 (blue), 200 (black) for a quadrupolar billiard (i?(0) = i?(l +ecos20)) of 
deformation e = 0.1. The eigenvalues are ordered with respect to their absolute value, 
and the horizontal axis denotes their relative order within all the eigenvalues (TV). 
Note that the unitary sector scales linearly with fc (corresponds approximately to the 
number of open classical channels which can be estimated to be 2 [fci?] ) . As the number 
of eigenvalues scales with the size of the system too, the unitarity border is identical in 
all cases, (b) Distribution of the eigenvalues in the complex plane for kR — 100. The 
solid line is the unit circle. In each case, the size of the system and hence the number 
of eigenvalues is iV = [6 x kR\ . 



For arbitrary shapes, K{s, s'; k) is however not unitary [31] and does not obey the 
spectral theorem. 

Nevertheless, a favorable property of this set is that for the finite-dimensional 
truncation of K{k), the spectrum can be roughly divided into a null space and a unitary 
sector (to be defined below). This can best be visualized by looking at the eigenvalue 
distribution of K. In figure[2]^a), we plot the absolute values of the eigenvalues {A*^*^(A;o)}- 
It's clearly seen that the distinction between null-space eigenvalues (|A''*''| ~ 0) and the 
unitary eigenvalues (lA*^*-*! ~ 1) becomes sharper for larger k i.e. in the semiclassical 
limit. At a typical value of k, the eigenvalues are distributed in the complex plane 
within the unit circle, and a fraction of the eigenvalues lie in the vicinity of the 
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unit circle representing the unitary sector (see figure [2](b)). The size of this unitary 
sector is approximately 2 [kR] , which corresponds approximately to the number of half 
wavelengths on the boundary [35] . 

Returning to the eigenvalue equation ([9]), we see that the quantization condition (IHl) 
can be rewritten as A(fcq) = 1. In other words, whenever we find an eigenvalue X^^\kq) 
at 1 + Oi in the complex plane, kg is a solution of ([H]) and u^'^^kg] s) is the associated 
quantized eigenvector. 

We will now argue that not only are the (unitary sector) eigenvectors of K{k) 
approximately orthogonal at a given k, but they also approximately diagonalize K{k) 
over a range 6kR ~ 0(1). (This range roughly corresponds to the one dimensional free 
spectral range of the billiard, which is roughly the number of wavelengths one can fit 
into the longest chord in the domain D. We will refer to it shortly as the "free spectral 
range".) Consider the eigenvectors calculated at two different but close values of the 
parameter k, say ko and ko + Sk. We can define the overlap between the eigenvectors 
calculated at these two different values by 

N 

{u^'\ko)\u^'\ko + 5k)) = J2{u^'^y{ko; si) u^^\ko + 5k; si). (11) 

1=1 

This operation is well-defined as long as we keep the system size constant. In figure [3l 
we start with an initial set of states \u^'^\ko)), i = 1, . . . , N and plot for subsequent 
k = [ko, ko + Ak] only the overlap of the various initial states with their maximal 
overlap partner. We would like to note that there is in general only one state at k 
that has a considerably larger overlap than all other states with an initial state u^*'*(/co)- 
Here we plot only a fraction of the initial eigenvectors for the sake of visibility, but this 
behavior holds in general over stretches 5kR ~ 0(1) of the parameter k for eigenvalues 
in the first and fourth quadrant of the complex plane (| argA*^*^| < 7r/2). The typical 
change in overlap over 5kR = 0.2 at kR ~ 100 is less than %1. 

An important consequence of this observation is that we can assign an identity to 
the eigenvectors even away from quantization [231 136]. To elucidate this point, consider 
the trace of one of the eigenvectors in figure HI The initial eigenvector is not quantized 
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(a) (b) (c) 

Figure 4. We follow an eigenfunction of the quadrupolar shape with e = 0.1 
via the method of highest eigenfunction overlap (see figure [3]). Starting from the 
quantized eigenvalue at kR ~ 20.725 wc follow it in steps of SkR = 0.025 through 
one quantization cycle up to kR — 24.175. (a) False color plots of the intensity of the 
traced wavefunction. (b) Eigenvalues at each snapshot. In red we trace the motion of 
the eigenvalue of the particular state plotted in (a) and in (c) we plot the corresponding 
error on the boundary given by = <fgjj ds\'tljfj,{r{s))\'^ . 



and we follow this state by the highest-overlap criterion over a range of 6kR ^ 3.5, 
a range that is larger than the free spectral range. We only plot here five instances 
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Figure 5. Eigenvalues of the quadrupole with e = 0.1. The speed d(j)/dkR of a number 
of eigenvalues around kR = 100. The corresponding eigenvectors are traced over one 
fourth spectral range, kR = [100.0, 100.2]. 



over which the state becomes quantized (top to bottom). An important feature of this 
behavior is the way an extra node is "pushed" into the bilhard. We have to note that 
over such large stretches of kR, an eigenvector typically undergoes avoided crossings. 
The avoided crossings happen predominantly around argA'-*-' ~ ivr in the second and 
third quadrants of the complex A-plane. This is the region of the complex eigenvalue 
plane where the null-space eigenvalues join the "unitary flow" (see figure [2](b) and the 
animation in the media section). However, the numerical method that we propose below 
utilizes the behavior in the first and fourth quadrants in the complex eigenvalue plane 
away from avoided crossings. 

A second key observation concerns the behavior of eigenvalues A^(/c) of K{k). 
This notation makes explicit the adiabatic identity of the eigenvectors that we have 
established above. With increasing k, the eigenvalue flow is counterclockwise. There is a 
clear distinction between the unitary eigenvalues which flow along the unit circle |A| = 1 
and the null-space vectors which accumulate at A ~ 0. The eigenvalues in transition 
that have an intermediate value of |A| follow a universal path (compare to the case of 
circular billiard in figure [TOj) and are added to the unitary flow at about (p = arg A ~ vr 
as noted above. This is the mechanism by which the density of states of the billiard 
eigenvalues increase, which according to the Weyl formula has the mean asymptotic 
behavior Pweyi{k) = kA/2n, where A is the area of the domain D. In figure [5] (a), we 
show that the phase speeds of the unitary eigenvectors, defined by v'^{k) = d(j)fj_{k)/dk 
is constant over a stretch of 6kR ~ 0(1). This is one of the main ingredients of the 
numerical diagonalization procedure that we propose in the next section. 
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Figure 6. Extrapolation error defined by equation (|14p for eigenfunctions of the 
stadium billiard of deformation L/R = 2 (see inset). The horizontal axis represents the 
initial value of the phase of the eigenvalue (/)° from which the solution is extrapolated. 
Data sets for kRo = 50, 100, 200, 300, 400 are plotted. 



4. An accelerated Fredholm root search and the accuracy of solutions 

Building on these observations, we propose the following numerical algorithm to compute 
both the billiard eigenvalues and the corresponding eigenfunctions (equation ([T])). We 
first determine the unitary eigenvectors, the eigenvalues and their corresponding phase 
speeds {Xfj,{ko),Ufj,{ko),v'^{kQ)} at a value k = ko equation ([9]). This requires two 
diagonalizations. We then extrapolate the quantization values k^ using the approximate 
constancy of the phase speeds 

K = ko + ^i^A2n-^l), (12) 

where (f)^^ = arglX^ko)]. The billiard eigenfunctions in the domain are then computed 
using the approximate u^ko) through 

Mr)=f dr'{s)G,{r,r\s)-k,)u,{k,-s). (13) 

JdD 

To assess the accuracy of the solutions we introduce the following quantity, the 
extrapolation error, 

p . ||(i-i^(MK(fco)ll2 .... 

\\Uf,{ko)\\2 

Here, || ■ II2 denotes the 2-norm. In figure El we plot the resulting error E^ko) for 
extrapolation from various values of initial fco- Instead of ko, we plot the error as a 
function of 0°. This provides a measure of the accuracy of the solutions as a function 
of the interval over which we extrapolate. This in turn determines the fraction of 
eigenvalues with a given accuracy. Note that a given </)|^ occurs at a different values of 
k for each fi. 

The data for different k^R in figure [6] demonstrates that despite the highly 
oscillatory nature of the higher lying excited billiard eigenfunctions, the error remains 
relatively constant as k is increased. A representative highly excited stadium state is 
plotted in figure [71 
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Figure 7. A representative plot of a quantized wavefunction for tlie stadium of 
deformation L/R = 2 at kR — 1300.02749 and an extrapolation error of 0.00296. 



We should note that the bilhard eigenf unctions presented here are domain- 
normahzed. As the normal derivative of the wavefunction u{s) on the boundary contains 
all information to determine the wavefunction throughout the domain, it's possible to 
express the normalization condition in terms of u{s) as [37] 

/ dsn-r{s)\u^iko;s)\^ = 2kl (15) 

JdD 

which then yields a ipfj,{x) which is normalized to unity in D. 

Next we compare the accuracy of the extrapolation method to that of SVD. In 
table [1], we compare the eigenvalues found via EVD extrapolation kj^^^ to those found 
by an SVD scan over an interval of [40.75, 50.25]. The extrapolation method can attain 
an accuracy obtained by an SVD scan at about 75 points, providing a factor of roughly 
10 in computation speed as seen in table [TJ Important to note is that the simple 
SVD scan will fail to account for all resonances. Only at a scan over 200 points have 
all resonances been resolved, which increases the factor to 270 The gain in speed at 

* As stated before, methods have been proposed to use several of the lowest singular values to resolve 
nearby resonances [20] , however we find that such an algorithm still leaves room for ambiguity at large 
wavevectors compared the the EVD extrapolation method. 
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Table 1. Comparison of the EVD and SVD method. The EVD extrapolation is 
performed at fco-R = 50. Column 1 contains the extrapolated billiard eigenvalues for 
a stadium of deformation L/R = 2. Columns 2-5 contain the eigenvalues obtained by 
an SVD scan at 25, 75, 200, and 5000 points in the interval [40.75,50.25]. Column 6 
contains the extrapolation error for the EVD eigenvalues in Column 1, and the final 
column contains the relative error of the EVD eigenvalues with respect to SVD5000. In 
the last row we have the average relative error for eigenvalues in Columns 1-5 compared 
to SVD5000. We also quote the computation time on a quad core CPU running at 
1.6GHz in the second row. Note that in this table, we show only a fraction of the 
eigenvalues computed. 



EVD 


SVD25 


SVD75 


SVD200 


SVD5000 






Rel. error 


10.41s 


35.22s 


102s 


272s 


6728s 










49.77668 


49.7700 


49.7700 


49.7725 


49.77180 


6.3505^ 


-3 


9.8116^;- 


5 


49.88509 


49.8900 


49.8900 


49.8850 


49.88560 


1.2706^ 


-2 


1.0193^- 


5 


49.89177 






49.8900 


49.88900 


3.7417S 


-2 


5.5620^- 


5 


49.94180 


49.9500 




49.9400 


49.93930 


1.6519S 


-2 


5.0195^- 


5 


49.94993 




49.9500 


49.9500 


49.94940 


7.4363S 


-3 


1.0725^- 


5 


50.03633 




50.0430 


50.0350 


50.03520 


1.1042S 


-2 


2.2719^- 


5 


50.04363 


50.0500 




50.0450 


50.04460 


5.5522S 


-3 


1.9338^- 


5 


50.08058 




50.0770 


50.0775 


50.07810 


1.0779S 


-2 


4.9640^;- 


5 


50.08856 


50.0900 


50.0900 


50.0900 


50.09020 


1.8785S 


-2 


3.2551^;- 


5 


50.16631 


50.1500 


50.1570 


50.1575 


50.15740 


5.5390S 


-2 


1.7773^;- 


4 


50.20171 


50.1900 


50.1970 


50.1950 


50.19410 


3.6940^ 


-2 


1.5175^- 


4 


50.23331 




50.2300 


50.2325 


50.23140 


7.6969^ 


-2 


3.8103^- 


5 


5.9723^-5 


1.7245^-4 


5.7066^-5 


1.1829£;-5 










a fixed accuracy will grow linearly with kR as the number of modes within a given 
interval of initial phases A0° will increase linearly with kR. We would like to emphasize 
that the accuracy of EVD method and the SVD method is in principle identical (this 
is clearly seen in comparing the minima attained in figures [9](a) and (b)) and it is the 
level of accuracy that is desired that will determine the speed enhancement obtained 
by the EVD method. We have implemented more complex extrapolation methods to 
provide a desired level of accuracy. The ultimate accuracy that can be attained scales 
exponentially with the number of discretization points on the boundary. This is shown 
in figure [H Finally, in table [21 we show the accuracy of the EVD method for a case 
where analytic solutions are available, namely the circular billiard. 

5. Relation to the SVD method 

In this section, we would like to clarify the relation between our method and the 
SVD method [201 ES]. In figure [H we compare the lowest few singular values cr^ik) 
to |1 — A^(/c)| which we find by diagonalizing K[k) at an arbitrary k within the spectral 
range plotted. We find that the plots are almost identical. This should not be surprising, 
because L{k) = 1 — K{k) is the matrix whose singular values are computed. A significant 
point is however that whereas the singular values cr^{k) are real (this is a choice of the 
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5 10 15 

rj 

Figure 8. Logarithmic plot of the error given by 1 1 — A^(fcp) | for a state of the stadium 
bilhard of deformation L/R ~ 2 quantized at kR « 20.2965 as a function of boundary 
discretization, rj is defined by iV = rikR 

Table 2. Accuracy of the extrapolation method for the circular billiard. We compare 
the solutions obtained by the EVD extrapolation method (A;^^^) to the solutions 
obtained analytically (fc^) by finding the j*'' zero of the Bessel function J,„(a;). We 
show the relative error (with respect to the analytic solution computed to a precision 
of 8) and the extrapolation error in the fifth and sixth columns, respectively. 



m 


j 




J.EVD 


Rel. error 




6 


1 


9.93610952 


9.93723495 


1.1327E-4 


1.8836580E 


-3 


1 


3 


10.17346813 


10.17983831 


6.2616E-4 


1.2763550^ 


-2 


34 


3 


49.95933191 


49.96813288 


1.7616E-4 


1.3162860E 


-2 


16 


9 


50.04460601 


50.05588673 


2.2541E-4 


2.1667890E 


-2 


85 


2 


99.98282066 


99.98522643 


2.4062E-5 


3.1056810E 


-3 


60 


8 


99.98510243 


99.98854423 


3.4423E-5 


6.0242240E 


-3 


24 


21 


99.99434362 


99.99440486 


6.1246E-7 


6.0912910E 


-4 


2 


46 


149.99854919 


149.98312449 


1.0283E-4 


3.0056170^ 


-2 


68 


19 


150.02814761 


150.01975466 


5.5943E-5 


1.4191580E 


-2 





48 


150.01188245 


150.02755229 


1.0446E-4 


3.2199340^ 


-2 


57 


23 


150.04477281 


150.04170015 


2.0478E-5 


4.9759550E 


-3 



numerical SVD routine) and obviously not analytic as a function of fc, \^{k) are complex 
(and can be shown to be analytic). These points can be put into a more formal setting 
by following the discussion in [38]. Considering additionally the singular vectors v of 
L\k) = 1 — K\k), it can be shown that a generalized complex singular value can be 
defined by a proper choice of relative phases of and f^, which is analytic as a function 
of k. These analytic complex singular values are exactly 1 — \^{k). The real singular 
value calculated by the numerical SVD routines are cr^(A;) = |1 — ^^l{.k)\. 
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Figure 9. (a) Comparison of the four lowest singular values crp(fc) of L{k) = 1 — K{k) 
and (b) the four eigenvalues A^(fc) of K{k) with the lowest |A^(fc) — 1|, for a range of 
kR = 27.25, . . . , 27.75 in a quadrupole with e — 0.3. 



6. Explicit results in the circular billiard 

In this section we will substantiate the above observations for an analytically 
solvable problem, namely the Dirichlet problem of a circular billiard. To solve the 
Laplace eigenvalue problem with Dirichlet boundary conditions for a circular quantum 
billiard analytically, we can we can write the Green's function, using Bessel addition 
theorems [391, as 



Go(r, r') = --H+(A;|r - r'\) = — Y^, H+ (A;r) J™(A;r')e^'"(*-<^ ) (16) 
for r > r'. Then, assuming r' is on the boundary and r is outside the circular domain, 

(17) 



Let us also Fourier expand the field 

m 

Evaluating the integral in ([3]) we are left with a diagonal kernel 

= {t7TkRRi{kR)3'JkR) - l) 5mm'- 

Thus, the singular values can be written as 

am{k) = 2- i7rkRR:^{kR)y^{kR) = 0. 
Using the Bessel identity 

J„(x)H^'(a;) - H^(a;)j;,(a;) 
this can be equivalently written as 

CTm{k) = J^(fci?)H+'(fci?). 



2i 

7TX ' 



(18) 

(19) 
(20) 

(21) 
(22) 



Note that the singularity condition yields the secular equation of the internal Dirichlet 
problem i.e. Jrn{kR) = and that of the external Neumann problem (kR) = (with 
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Figure 10. The red curve shows the parametric behavior of the complex eigenvalues 
of the circle Am(fc), via equation ((23|) with m = 10 and kR = 0..100. In black circles 
we show the eigenvalues found numerically for kR — 100 with iV = [10 x kR]. 



Sommerfeld radiation conditions). The latter doesn't have any solutions on the real axis 
while the former has all its solutions strictly on the real axis. 

Now let us look at the eigenvalue problem and the extrapolation method. The 
eigenvalues are parametrically dependent on k and given by 

Xm{k) = -1 + inkR}l^{kR)J'^{kR). (23) 

In figure [10] we show this parametric behavior reproduces the general features observed 
for smoothly deformed shapes (compare to figureE]), in particular the transition behavior 
of eigenvalues from the null-space to the unitary sector. 

Using Debye asymptotic expansions of the Bessel Functions, one can show that for 
m < kR (m ^ 1) 

Xm{k) ~ e**, (24) 

where $ = 2m(tan/5-/3)+7r/2 and cos^ = m/x. For m > kR, \Xm\ ~ q-'^HP-^^^ P) <^ i. 
Note that the transition region around m ~ kR (which corresponds to the behavior in 
the transition region) is not represented uniformly by the above expressions. We thus 
find that the speed of the unitary eigenvalues (in this case m < kR) are asymptotically 
given by 

v'l'ik) ^2sm(3. (25) 

Hence, the change in speed is asymptotically small in kR {dv^/dk ~ (kR)^^) as is 
observed numerically for arbitrary smoothly deformed shapes. 

7. Conclusion 

We have presented an efficient and robust algorithm to calculate eigenvalues of the 
Laplace operator based on a novel Fredholm formulation. We have shown that 
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approximately of the order of kR eigenvalues can be found with just two diagonalizations 
and no root search. This overcomes a formidable problem faced by diagonalization 
algorithms based on SVD for finding large eigenvalues: distinguishing real from false 
minima in singular values when the density of states pweyiik) is large. 
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